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ABSTRACT 


GRID GENERATION AND FLOW 
COMPUTATION ABOUT A MARTIAN 
ENTRY VEHICLE 


A number of vehicles are currently being proposed for a manned mission 
to Mars. One of these vehicles has a modified blunt-nosed cone configuration. Ex- 
perimental results have been obtained for this vehicle in 1968. These results show 
lift-over- drag ratios comparable to those needed for Mars entry. Computations are 
performed here to verify the earlier results and to further describe the flight charac- 
teristics of this vehicle. 

An analytical method is used to define the surface of this vehicle. A single- 
block volume grid is generated around the vehicle using an algebraic Two- Boundary 
Grid Generation algorithm (TBGG) and transfinite interpolation. Euler solutions 
are then obtained from a Langley Aerothermodynamic Upwind Relaxation Algorithm 
(LAURA) at Mach 6.0 and angles of attack of 0, 6, and 12 degrees. 

The lift coefficient determined from the LAURA code agrees very well with 
the experimental results obtained in the 1968 study of this vehicle. The drag and 
pitching moment coefficients, however, are underestimated by the code since viscous 
effects are not considered. Contour plots of the flowfield show no evidence of separa- 
tion for angles of attack up to 12 degrees. 
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Chapter 1 

INTRODUCTION 


Aviation and aerodynamics have developed rapidly since the Wright brothers 
first flight in 1903. In the past decade, man has gone out into space and returned to a 
conventional landing on earth (via the space shuttle), spent long durations in space, 
and sent unmanned space probes beyond our solar system. The next step in the 
evolution of space travel seems to be the manned exploration of the planets. Mars, 
because of its close vicinity to earth and its physical characteristics, appears to be the 
first choice for such a mission. This prompted the study of the atmosphere of Mars in 
order to determine what type of vehicle and mission could be used to land there. One 
possibility seems to be a mission in which a vehicle enters the atmosphere, slows to 
low supersonic speeds, and deploys a parachute. The vehicle can then use its engines 
to slow down, release the parachute, and land on Mars under its own power. 

Many of these types of vehicles were studied in the 1950’s and 60’s for use 
in missions to the moon. These vehicles are now being reexamined to see if they may 
be used in a mission to Mars. One of the vehicles being examined comes from a 1968 
technical paper entitled “Aerodynamic Characteristics of a Modified Cone-Conical- 
Frustum Entry Configuration at Mach 6.0” [1]*. This paper contains a description of 
the vehicle and results of experimental tests performed on the vehicle. A review of 
these results shows that this vehicle might be useful in a manned mission to Mars. 

’ The numbers in brackets indicate references. 
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Lift-over-drag ratios are comparable to those needed for entry into the Mars atmo- 
sphere at hypersonic and then supersonic speeds. 

Methods used to determine the flow characteristics of aircraft and spacecraft 
have changed a great deal since 1968. Many of the advances made in the areas of avi- 
ation and aerodynamics were made possible by the advances made in Computational 
Fluid Dynamics (CFD). CFD now plays a dominant role in the aerospace field because 
of its effectiveness as a design tool and as a compliment to experimental tests. The 
last two decades have brought the introduction of flow solvers which are capable of 
solving the partial differential equations of fluid motion quickly and efficiently. These 
solvers have been verified experimentally and are now widely used in the design of 
aircraft and spacecraft. Continuing improvements in high speed and large memory 
digital computers also act as a catalyst for the use of CFD in the future. The purpose 
of this study is to verify the 1968 results and to gather new information about the 
aerodynamic characteristics of the vehicle. 

A typical CFD problem is divided into two steps - grid generation and flow- 
solution. The surface grid in this study is constructed analytically from the original 
model details [1]. These details, however, lack some information needed to completely 
define the nose of the vehicle. Therefore, a smooth surface grid is constructed by 
making some assumptions at the nose and following the original model details as 
closely as possible. 

The volume grid can be constructed using many different methods [2,3]. 
Some of these include transfinite interpolation, multi-surface interpolation, elliptic 
grid generation, and hyperbolic grid generation. Transfinite interpolation is used in 
this study because of the control it allows over the volume grid, its computational 
efficiency, and the ease in which it can be implemented. Transfinite interpolation 
requires that the six sides of the volume grid be created before interpolation be- 
gins. Two of the six sides are created analytically while the remaining four sides are 
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generated using an interactive algebraic two-dimensional grid generation algorithm 
called Two-Boundary Grid Generation (TBGG) [4]. Interpolation is then performed 
to compute the interior of the volume grid. 

Flow solutions can be obtained to varying degrees of accuracy by many 
different methods. Finite-difference, finite- volume, and finite-element techniques are 
widely used as a means of solving the partial differential/integral equations of fluid 
motion. These equations in their viscous, compressible form are the Navier-Stokes 
equations. Flow solutions obtained from these equations are often computationally 
expensive. They can, however, be modified in such a way that the diffusion terms 
are discarded [5]. These modified equations are called the Euler equations. The 
solution of these equations is computationally more efficient than the solution of 
the full Navier-Stokes equations and will often produce good estimates of lift and 
pressure distributions for the vehicle. Values of drag and pitching moment, however, 
are underestimated, since only form drag and not viscous drag is included. 

In this study, the Langley Aerothermodynamic Upwind Relaxation Algo- 
rithm (LAURA) [6,7,8] is used to solve for the flowfield around the vehicle. LAURA 
is a viscous code which has been modified to solve the Euler equations. A literature 
survey gives a comparison of numerical and experimental results [9] validating the 
code, and there are a number of examples of the use of the LAURA code for different 
configurations [10,11]. 

The experimental tests performed on the vehicle in 1968 [1] were conducted 
at Mach 6.0. In the present study, computational tests are also performed at Mach 
6.0 in order to verify these results and to further describe the flow around the vehicle. 
A free-stream Mach number of 6.0 falls within the hypersonic range of fluid flow. 
The characteristics of blunt-nosed bodies in hypersonic flow is well documented [12]. 
The flow is characterized by the formation of a strong shock wave, subsonic flow on 
the leading edge of the blunt nose behind the shock, and supersonic or hypersonic 
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flow behind the shock in the surrounding regions. A compression, overexpansion, and 
recompression is also expected near the nose in the streamwise direction. Similarly, 
hypersonic flow around conical bodies is also well documented [12]. It is characterized 
by strong shock waves coming off the cone at an angle which can be determined 
analytically based on cone angle and free-stream velocity. 

This study is divided into chapters which follow a logical sequence. A de- 
tailed description of the generation of the surface and volume grid is given in Chap. 2. 
The inviscid governing equations for compressible fluid flow are given in Chap. 3. The 
LAURA code and its modifications to compute inviscid flowfields are also described 
in this chapter. Chapter 4 gives a summary of the results and comparisons between 
experimental and computational results. The concluding remarks and the future 
directions of this study are presented in Chap. 5. 


Chapter 2 

GRID GENERATION 


2.1 Introduction 

The configuration of the spacecraft is taken from a 1968 NASA Technical 
Note [1], This design evolved from a basic body of revolution of a 15.07° conical 
forebody, 11.3° conical afterbody, and a spherical nose. Figures 2.1 and 2.2 show the 
details of the original model [1]. The afterbody consists of two cylinders mounted in 
parallel to the exterior of the 1 1.3° inner cone. These cylinders are tangentially faired 
into the inner cone as shown in Fig. 2.1. The forebody consists of two 6.64° cones 
mounted externally to a 15.07° inner cone. These outer cones are also tangentially 
faired into the inner cone. The spherical nose then fits on the forebody such that 
C l (first derivative) continuity is maintained. Figure 2.3 defines the location of the 
nose, forebody, and afterbody. 


2.2 Surface Definition 

The surface grid is created analytically from the original model details [1]. 
Lines of constant £ are constructed as lines of constant angle 9 in the physical space 
where 

9 = arctan — . (2.1) 

x 
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Fig. 2.1 Original model details (all dimensions in cm.). 
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Fig. 2.3 Nose, forebody, and afterbody locations. 
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Lines of constant r/ are constructed as lines of constant axial distance z in the physical 
space. Therefore, <f-lines lie in the x-y plane and outline cross-sections of the model 
in a plane perpendicular to the model’s axis of symmetry. It is necessary to find the 
outline of each of these cross-sections in Cartesian coordinates to define the surface 
of the model. 

The afterbody contains cylinders mounted parallel to the walls of an 11.3° 
inner cone. The cross-sectional cut of a cone perpendicular to its axis of symmetry 
creates a circle in the x-y plane (Fig. 2.4a). A cross-sectional cut of a cylinder at 
11.3° creates an ellipse (Fig. 2.4c). The minor axis, a„, of this ellipse is equal to the 
radius of the cylinder. The major axis, b a , is equal to 


t r _fvj 

a cos 11.3° 


( 2 . 2 ) 


The forebody contains two 6.64° cones mounted at an angle of 15.07° to 
the axis of the inner cone. A cross-sectional cut of these cones at an angle of 15.07° 
creates an “egg-shaped” ellipse (Fig. 2.4b). The length of the major axis on one side 
of the ellipse is different than that of the other. Noting that the radius of the inner 
cone is always larger than that of the outer cone for each cross-section (Fig. 2.1), the 
innermost major axis is neglected. The tangent between the inner cone and outer 
cone will always intersect the outermost side of the egg-shaped ellipse. The outer 
major axis is calculated as 


b } = 


(2.3) 


cos 15.07° 

The radius r oc is the radius of the outer cone at the point of intersection of the outer 
cone axis and the cross-sectional plane (Fig. 2.4b). The minor axis, a/, is equal to the 
radius of the outer cone at this point. The spherical nose begins at the point where 
the outer cone reaches its apex. A cross-section of the spherical nose creates a circle 
in the x-y plane. 

Figure 2.5 represents a typical cross-section of a region in the forebody or 
afterbody. The radius, major axis, minor axis, and ellipse center are now determined 
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H 



Fig. 2.5 Forebody or afterbody cross-section. 
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from the information given above. This information is necessary in order to determine 
the points of tangency on the circle and ellipse. These points, along with the first 
derivatives at these points, represent six unknowns. Six equations can be constructed 
to determine these unknowns. The equations are: 

Equation of a circle 


X c 2 + Vc 2 = 


(2.4) 


Equation of an ellipse 
Equation of a triangle 


(2.5) 


r, c 2 + [(x e - 1 ,) 2 + (y e - y.) 2 | = [x 2 + y] 


( 2 . 6 ) 


Equation of a tangent at (x c ,y c ) 


dye _ _£c 
dx c y c 

Equation of a tangent at (x e ,y e ) 

9y e ( Xe 3-ec \ ^ 

dx e ~ V ye ) P 


(2.7) 


( 2 . 8 ) 


Equal tangents equation 

dye _ dy^ 
dx c dx e 

Equations (2.4 - 2.9) can be simplified into two equations 


f(x c , x e ) — 




a 


fl - [(Xe - X ec )/6] 2 * 2 




= 0 


( 2 . 10 ) 


J (x I ,x«)=2r„>-2x e x.-2,.{[l-(^ T ^) J ][r, t J -x c 2 ]}'/ 2 = 0 , (2.11) 

where b is the major axis of the ellipse at each cross-section (b a or b/), and a is the 
minor axis of the ellipse at each cross-section (a a or a/). Equations (2.10) and (2.11) 
can be combined into one non-linear equation. It is more convenient, however, to 
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leave them as separate equations and solve for x c and x e using a Newton- Raphson 
method. The matrix form of this method is 


» 

h' 

o 

[ Ax c 


1 

'-s 

L 9x c 9x e J 

%> 

H 

<1 
1 


1 

i 

< 


where Ax c =x c n+1 — x c n , 

Ax«.=x e n+l - x e n , 


( 2 . 12 ) 


and subscripts x c and x e denote differentiation with respect to these variables. Initial 
values x c n and x e n are required to begin the solution of Eq. (2.12). These values can 
be found by simply estimating the points of tangency on the circle and ellipse. These 
estimations, however, can be inaccurate and in some cases result in the divergence 
of the solution. To avoid this problem, a new method is developed to find the initial 
values. As noted earlier, the point of tangency on the ellipse is always on the side 
away from the inner cone. In other words, x e is always between x ec and x ec +b. This 
allows an initial value of x ec + .5b to be used for x e ". The initial value x c n is then 
found by solving Eq. (2.10) for x c as a function of x e , i.e. 


(x e 

- x ec ) 2 a 2 r tc 2 / I 

( *1 ( X r Xec 

) k \ 


b ) 


u- 

(x e - x ec ) 2 a 2 / 



iJ 


(2.13) 


The Newton- Raphson method works well in finding x c and x e over most of 
the model. Near the nose, however, the solution fails to converge. The bisection 
method is then used to determine x c and x e over the remaining cross-sections. This 
method requires that there be only one unknown. Equation (2.13) is therefore used 
to find x c as a function of x e . This equation is then substituted into Eq. (2.11) such 


that 

g = g(x e ) . ( 2 - 14 ) 


The solution is initially bracketed by x ec and x ec + b. The value of x e determined 
from this method is then substituted into Eq. (2.13) to determine x c . 
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The coordinate lines £ and r\ can now be constructed on the surface of the 
model in the physical space. The y-coordinates corresponding to x c and x e are found 
from Eqs. (2.4) and (2.5), respectively. These equations, along with the equation 
of a line between (x c ,y c ) and (x e ,y e ), are used to create the coordinate line £ in the 
physical space. Coordinate line rj is constructed using Eq. (2.1), noting that each tj 
line is at a constant angle 0 in the physical space. 

A side and top view of the surface grid is given in Fig. 2.6. These views show 
concentrations of grid lines at certain locations on the model. The concentrations are 
necessary in order to completely resolve possible shocks in these regions. This will be 
discussed in the following chapter. If X is a uniformly spaced coordinate, the equation 


Y = 


„KX 


- 1 


t K - 1 


(2.15) 


where 


0< Y <1, 
0< X <1, 


is used to concentrate grid points near Y=0 for K> 0, or near Y=1 for K< 0. High 
concentrations of grid points result when large positive or negative values of K are 
used. Successive concentrations are obtained by connecting a number of curves cre- 
ated by Eq. (2.15) such that they are C 1 continuous. 


2.3 Nose Discontinuity 

The construction of the afterbody is fairly straightforward. Figure 2.1 gives 
all the information necessary for its construction. The information on the forebody 
and nose, however, is incomplete. Specifically, more information is needed to de- 
termine how the spherical nose fits onto the forebody. Since this information is 

unavailable, some assumptions are made. 

The initial assumption is that each outer cone reaches its apex at the same 
axial location that the inner cone is truncated (Fig. 2.7). A cross-section of the 
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Fig. 2.6 Side and top view of surface grid. 










model taken perpendicular to the model’s axis of symmetry at this point will result 
in a circle. Theoretically, a spherical cap can then fit onto the model. Upon close 
examination of the model, however, this method results in a spherical nose that has a 
point discontinuity in slope at the apex of each outer cone. The outer cones approach 
the spherical cap at an angle of 25.59° to the model’s axis of symmetry while the 
inner cone approaches at an angle of 15.07° (Fig. 2.7). If the spherical cap is made C 1 
continuous on the inner cone, then it will not be C 1 continuous at each outer cone’s 
apex and vice versa. 

A logical choice is to make the spherical cap C 1 continuous on the inner cone 
and to smooth out the two point discontinuities in slope caused by the outer cones. A 
smoothing algorithm is used to smooth these points locally. Details of this algorithm 
can be found in Refs. 11 and 13. The algorithm creates piecewise continuous cubic 
splines through discontinuous regions allowing the user to smooth particular regions of 
the model without affecting the surrounding regions. When applied to the model, the 
program essentially “fills in” the step change in slope at each discontinuity (Fig. 2.8). 
This, however, creates a new problem. The regions adjacent to the smoothed region 
are now “lower” than the smoothed region as a result of the program filling in the 
discontinuity. This causes “valleys” to appear in the surface. These valleys cannot 
remain on the finished model. They can essentially be flattened out by smoothing 
a larger region of the model. This, however, is not done since the original model 
configuration is lost when large regions are smoothed. If a reasonable model is to be 
created, it will be necessary to reexamine some of the original assumptions made at 
the nose. 

The assumption that the outer cones reach their apex at the same axial 
location that the inner cone is truncated, is modified slightly to allow a nose to be 
created that does not require smoothing. This is done by keeping the spherical cap 
C 1 continuous on the inner cone, but allowing each outer cone’s apex to intersect 
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Fig. 2.8 Filling in discontinuity. 
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the spherical cap such that they are C 1 continuous with the cap (Fig. 2.9). This 
essentially moves each outer cone towards the spherical cap slightly, allowing their 
apexes to intersect the cap tangentially. In the small region where the outer cones 
extend beyond the forebody, the outer cones are tangentially faired into the spherical 
nose rather than the inner cone. This requires that the inner cone radius, r tc , be 
replaced by the spherical nose radius, r jn , in Eqs. (2.4) and (2.6), and subsequently 
in Eqs. (2.10) and (2.11). These equations are then solved using the same methods 
described in Sec. 2.2. 

The changes imposed on the nose of the model have an effect on the forebody- 
afterbody junction of the spacecraft. The 11.3° inner cone of the afterbody and the 
15.07° inner cone of the forebody are still C° continuous. No changes are made on 
either inner cone when the nose is changed. The outer cones of the forebody and 
the outer cylinders of the afterbody, however, are no longer C° continuous. If the 
locations of the apex of each outer cone is changed, the location of the entire cone is 
subsequently changed. The C° continuity is maintained at the junction if the 6.64° 
outer cones are replaced with cones of frustum angles less than 6.64°. This, however, 
is not done since the outer cones of the forebody are explicitly given to be 6.64° in 
the original model details (Fig. 2.1). Instead, the angle of inclination of the outer 
cones is changed to accommodate the changes made in the nose. 

Figure 2.10 shows the new configuration of the forebody- afterbody junction. 
The difference in radius between the base of the outer cones and the outer cylinders 
causes a small C° discontinuity at the junction but this can easily be smoothed using 
the smoothing algorithm described earlier [13]. Figure 2.11 shows the discontinuity 
before and after smoothing. Figure 2.12 gives a detailed schematic of the new model 
details and the two programs which generate the surface grid are given in Appendix A. 









Fig. 2.11 Surface discontinuity: (a)before smoothing, and (b)after smoothing 
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2.4 Volume Grid 

A single-block volume grid can now be generated about the surface grid. 
The volume grid is a C-0 type grid with coordinate direction £ defined normal to 
both £ and rj in the computational domain. Flow solutions for the vehicle are obtained 
at various angles of attack but at 0° yaw. It is therefore only necessary to generate a 
volume grid around half of the vehicle. The vehicle enters the atmosphere with the 
x-axis (Fig. 2.12) oriented horizontally. The vehicle is therefore divided along the 
y-z-axis and a grid is generated around one half (side view of Fig 2.6). 

The volume grid is generated in two steps. First, the top symmetry, bottom 
symmetry, side, exit, and outer grids are created. Then the interior of the volume grid 
is generated using transfinite interpolation. The top symmetry, bottom symmetry, 
side, and exit grids are created using an interactive, algebraic, two-dimensional grid 
generation program called TBGG (Two-Boundary Grid Generation). This program 
is described by Smith and Wiese in a 1986 NASA technical paper [4). TBGG first 
reads an input data file containing points which define the four outer boundaries of 
the 2-D grid. The points on the first boundary of each 2-D grid are defined by the 
existing surface grid. These points are fixed in TBGG so that they always coincide 
with the points on the surface grid. The three remaining boundaries of the 2-D grids 
are determined in such a way as to capture the shocks coming off the vehicle. 

The top, bottom, and side grids share one common boundary. This bound- 
ary is along the negative z-axis and is shared by all 2-D grids in the r)-£ computational 
plane. This boundary is a polar singularity and lies along the stagnation line of the 
vehicle at 0° angle of attack. The length of this line should be at least 1.5 times 
the shock stand-off distance expected in front of the blunt nose in order to properly 
capture the shock. At Mach 6.0 the shock stand-off distance is approximately .006 d 
( d being the base diameter of the 11.3° inner cone shown in Fig. 2.12) for a .08 d 
diameter sphere [14]. It is very difficult to create a reasonable grid within a .009 d 
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region between the blunt nose and the volume grid’s outer boundary. The grid lines 
coming off the surface would require a high degree of curvature in order to maintain 
orthogonality at the surface. Therefore, the outer boundary is moved away from the 
blunt nose a distance of three times the diameter of the spherical nose. The increased 
distance allows a more uniform grid to be created in front of the nose. 

The two remaining boundaries of the top, bottom, and side grids are de- 
termined from the shock wave angle created by a cone. The shock wave angle for a 
15.07° cone at Mach 6.0 is approximately 19° while the shock angle for a 25.59° cone 
is approximately 30° [15]. Therefore, in order to capture the shock wave around the 
vehicle, the outer boundaries of the top and bottom grids are inclined to an angle 
of 25° while the outer boundary of the side grid is inclined to an angle of 35°. If 
computations are performed at different angles of attack, the outer boundary of the 
top grid must be inclined to an angle greater than 25° in order to properly capture the 
shock. Estimating that the shock wave angle for a 15.07° cone at an angle of attack 
of 15° is equal to the sum of the shock wave angle (at 0° angle of attack) and the 
angle of attack, an outer boundary inclined at 40° will capture the shock. Therefore, 
this outer boundary captures all shocks for angles of attack up to 15°. Using this 
information, the four boundaries of the top, bottom, and side grids can be generated 
for input into TBGG. 

The outer boundaries of the top, bottom, and side grids are used in the 
creation of the boundaries of the exit grid. The top, bottom, and surface grids 
coincide with the exit grid on three of the four boundaries. The point distributions on 
these boundaries are fixed in TBGG so that they will always coincide. The remaining 
boundary is defined by the intersection of the top, bottom, and side grid with the 
exit grid. This boundary is essentially a “stretched” circle. The radius of the circle 
is determined by linear interpolation between intersections of each grid with the exit 
grid. 
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TBGG begins with the outer boundaries of each 2-D grid. The user then 
specifies the number of grid points, point distributions (assuming they are not already 
defined as described above), orthogonality, etc. in order to generate a 2-D grid. The 
grid is then viewed and modified until a satisfactory grid is created. Point concen- 
trations are made at the nose, near the vehicle, and near the forebody/afterbody 
junction in order to resolve high gradients expected in these regions. The resulting 
2-D top, bottom, side, and exit grids are shown in Fig. 2.13. 

The next step in generating the volume grid is to create the sixth side of 
this grid. The polar singularity, surface, top, bottom, and exit grids are the five sides 
of the volume grid which now exist. The side grid is not a side of the volume grid 
but an intermediate grid between the top and bottom grids. It is used to control 
the shape and point distributions of the volume grid in the physical space. The final 
side of the volume grid is called the “cap”. The cap is a 3-D surface in the physical 
space but exists as a 2-D surface in the computational space. This grid is in the £-77 
computational plane opposite the surface grid. Two boundaries of the cap are defined 
by their intersection with the top and bottom grids. The third boundary of the cap is 
defined by the location of the polar singularity while the fourth boundary is defined 
by the exit grid. All of these boundaries have point distributions which are defined in 
the physical space (excluding the polar singularity which is a point). It is necessary 
to define the interior of the cap before generating the interior of the volume grid. 

TBGG is not used to generate the interior grid of the cap since the cap is a 
surface in the physical space. Instead, the interior of the cap is generated based on the 
point distributions given on it’s boundaries and on the point distributions given on the 
side grid. The side grid intersects the cap along an //-coordinate line lying between the 
top and bottom grid boundaries. The points on this line are used to control the shape 
of the cap which in turn controls the shape of the volume grid. The cap is generated by 
first considering cross-sections of this grid projected on the x-y plane. Each coordinate 
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Fig. 2.13 Top, bottom, side, and exit grids. 
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line £ of the cap is generated in the same way that the outer boundary of the exit 
grid is generated. Three points are defined on this line by intersections of the cap 
with the top, bottom, and side grids. Linear interpolation between these points is 
used to find the radius of a stretched circle in the x-y plane. The circle is generated 
through these points to create each ^-coordinate line. The point distributions along 
the £ coordinate lines are found based on the angular distributions of the points along 
the joint boundary of the cap and the exit grid. Therefore, each coordinate line y lies 
at a constant angle 6 , where 8 is defined by Eq. (2.1). Figure 2.14a shows the shape 
of a typical cap cross-section projected on the x-y plane. 

The x and y locations of the cap’s interior are now defined in the physical 
space. The z locations of the grid are found by looking at the grid projected on the 
y-z physical plane (Fig. 2.14b). Three points are defined on the £ coordinate line by 

f 

the intersections of the top, bottom, and side grids with the cap. A C 1 continuous 
line must be constructed through these points in order to define the z locations of 
the ^-coordinate lines in the physical space. This is done by creating a quadratic 
equation of the form 

z = a x y 2 + a 2 y + a 3 , (2.16) 

where aj, < 12 , and a 3 are coefficients found by solving Eq. (2.16) simultaneously 
through the three points. This equation, therefore, creates a quadratic curve which 
is C 1 continuous through these points. Since the y locations of each coordinate line 
77 are known, the z locations of these coordinate lines are also found from Eq. (2.16). 
The interior of the cap is therefore completely defined in the physical space and is 
shown in Fig. 2.15. 

Transfinite interpolation is now used to generate the volume grid. The 
theory of transfinite interpolation as described by Gordon and Hall [16] is a general 
concept of multivariate interpolation between any number of surfaces. The single- 
block volume grid needed for this vehicle requires interpolation between the two 





Fig. 2.15 Definition of cap in the physical space 
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outermost planes of the £, t?, and £ coordinate directions. The transformation function 
defines the transformation from the unit cube computational domain to an 
arbitrarily-shaped region in the physical domain. Using this definition, the transfinite 
interpolation procedure is a recursive algorithm which is used to generate the interior 
of the volume grid. This algorithm is given as 

/(£,»?, C) = M£,^C),y(f, 7 ?,C)iz(£,77iC)] T i 

fx(t,v,0= 0 , 

h=\ 

= hit i*7iC) + yjM MU) ~ fiiZiVhiQ] ' 

h=i 

f((,rhO= [chit,!) ~ Mt,V^h)} • ( 21T ) 

h=l 

The planes a, 6, and c are defined as 

ai(»7,C) = /(0» T ?>C) > a l(j;,C) = /(U'bC) 1 

Su.o-/Uo,o , fi«,c) = /teu) , 

*(*,»?) =/tf.M) > <S(£,7) = /U>M) • ( 2 - 18 ) 

The univariate blending functions a, /?, and 7 are of the form 

ai(O = 1- 0(O » “2(0 = 0(0 ’ 

Pi(rj) = 1 - ^(17) , 02(T)) = H T l) . 

7.(0 = 1 - WC) . 7.(0 = WO , < 2 - 19 ) 

where t/>(X) = , 0 < X < 1 

The exponential equation above is identical to Eq. ( 2 . 15 ) of Sec. 2 . 2 . The effects of 
different values of K on are described in that section. 

Transfinite interpolation is performed on the volume grid in two steps. First, 
interpolation is performed between the top and side grids with the top grid defined 
as d[ and the side grid defined as 02. Then, interpolation is performed between the 
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side and bottom grids with the bottom grid defined as <xi and the side grid defined as 
02 - C l continuity cannot be guaranteed across the side grid boundary if the volume 
grid is generated in this manner but visual inspection has shown that a high degree 
of continuity is maintained across this boundary. 

A value of 2 is chosen for K in Eq. (2.19). Noting the definitions of a, /?, 
and 7 in Eq. (2.19), and the definition of the transformation function in Eq. (2.17), 
the majority of the emphasis is placed on the a j, b\, and c\ surfaces when K 2. This 
causes the point distributions on the surface, nose, top, and bottom grids to have the 
major influence over the distribution of grid points within the interior of the volume 
grid. A few intermediate grids are shown in Fig. 2.16 to demonstrate the results of 
the transfinite interpolation. 

The grid is constructed with 101 points in the ^-direction, 51 points in the 
^-direction, and 35 points in the ^-direction. This grid size is suitable for Euler 
calculations which are discussed in the following chapters. Figure 2.17 demonstrates 
the mapping from the physical domain to the unit cube computational domain. 
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Fig. 2 .16 Intermediate surfaces resulting from transfifiite interpolation 



Chapter 3 


GOVERNING EQUATIONS AND 
SOLUTION TECHNIQUE 

3.1 Introduction 

There are a number of methods currently being used to determine the flow 
characteristics around aircraft and spacecraft. The method used in this study is based 
on an upwind-biased finite volume algorithm developed by Gnoffo [6,8]. The viscous 
code based on this algorithm is called the Langley Aerothermodynamic Upwind Re- 
laxation Algorithm (LAURA). This code has been modified by Weilmuenster. Smith 
and Greene to compute inviscid flowfields [11]. A brief description of the LAURA 
code and the implementation of the inviscid boundary conditions are given here. 

3.2 Governing Equations 

The governing equations for this code are the three-dimensional Euler equa- 
tions for a compressible perfect gas. The integral form of these equations is [11] 

J J J q t d.0. + J J h da = 0 . (31) 

Expressing Eq. (3.1) in finite- volume form for a single six-sided cell (Fig. 3.1) in the 
computational domain gives 


35 




37 


.j,* = - 






+[(E^)ij+i k ~ + [( Gcr )i,j,*:+i ~ ( G<7 ).,j,*:-i]} (3-‘~) 


where ; St — ^n+i — t n 

The variable fi is defined here as the cell volume while cr is the cell face area. The 


vector q is defined as 


where 


q = 


■ 

p 

pu 

pv 

pw 

pE t 


E t = e + 1/2 (u 2 + v 2 + w 2 ) 


(3.3) 


(3.4) 


The variables p and e in Eqs. (3.3) and (3.4) are the nondimensional values of density 
and internal energy per unit mass, respectively. The variables u, v, and w are the 
Cartesian components of velocity in nondimensional form. The subscripted lower-case 
letters indicate cell center values, unless offset by a half, in which case they indicate 
values at a cell face (Fig. 3.1). 

The normal inviscid fluxes at a cell face (E, F, G) all have the form of H 


shown below 


H /+ i = - [a I (qi)’ + fri (qi+i)" 

■ ( 3 - 5 ) 

I/ ( + l 2*2 2 

The asterisk represents the time level and may be either the n or n-f-1 level, depending 
upon the cell face being evaluated. If variables with this superscript happen to contain 
information referenced in the i, j, k cell center, these values are linearized using 


K 


n+1 __ 


_ dK 

K +ir 

dq 




(3.6) 
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where K is a dummy variable representing the value to be linearized. If these lin- 
earizations are not included, an explicit rather than an implicit algorithm is created. 
The functions a w and b w are weighting parameters defined in terms of cell volumes as 


= 2 


n 


/+ 1 


fb + fb+i 

h = Q gj 

“ Vi + n i+1 


(3.7) 


The inclusion of the a w and b w parameters lessens the effects of grid stretching near 
the axis singularity and in other regions where the grid is highly stretched [7]. The 
choice of a w and b w is empirical, and other formulations are suggested in Ref. 6. 
For any cell face, the inviscid normal flux, I, is computed using Roe’s averages of 
cell-centered values and has the form 


1 = 


P U 

pull + Pn x 
pvU + Pn y 
pwU + Pn z 
(pE, + P)U J 


(3.8) 


In this equation, P is the nondimensional pressure and U is the contravarient velocity 
normal to the cell face. Variables such as the unit vector normal to a cell face, 
n = n x i + n y j + n t k, inverse distance between cell centers, v, eigenvalue matrix, 
A, and the right and left eigenvector matrices, M~* and M, can be found in Ref. 6. 
This reference also contains the definitions of the unit vectors tangent to a cell face, 
the cell volume, face area, and timestep. 

A first-order-accurate flux is computed when 6 is equal to zero in Eq. (3.5), 
while a second-order-accurate flux is computed when 9 is equal to one. Because the 
two terms are explained in detail in Ref. 6, only a brief outline is given here. The 
flux shown in Eq. (3.5) can be thought of as a second-order approximation of the 
flux at a cell face (first two terms), minus a dissipation term (remaining terms). If 
this dissipation is not included, the algorithm is equivalent to a centrally-differenced 
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algorithm. The first-order dissipation is given as 

-M |A| s* , 
v 

where 


s* = M"'Aq* 


(3.9) 


The change in q across a cell face, Aq, is computed via the upwind differencing 
method attributed to Roe [17]. The variable u is included here for the same reason 
a and b are included in Eq. (3.5). The matrix A contains the absolute values of the 
eigenvalues. Roe’s first-order dissipation term is the exact solution to the approximate 
initial value Riemann problem which is 


q t + JAq = 0 . (3.10) 

These three matrices (M, M - ', and A) relate to the inviscid flux vector Jacobian J 
by the equation 


J = = MAM - ' 


(3.11) 


Roe’s first-order dissipation term can be thought of as two flux differences across a 
cell face - one flux difference associated with the positive eigenvalues and the other 
associated with the negative eigenvalues, such that 


M |A| M“' Aq* = AI* - AI" = |J| Aq* . (3.12) 

Equation (3.12) is an exact equality if the matrices are evaluated using the Roe’s 
averaged variables. For a nondimensional case, Roe’s method interprets the inviscid 
zone of dependence correctly. Although not provided, the same is assumed for a 
three-dimensional flow. Allowing this assumption, if a flow is supersonic and has all 
positive eigenvalues, a flux is constructed based entirely on the upstream information. 
Problems with eigenvalues of mixed signs are handled accordingly. 
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Second-order accuracy is achieved with the Symmetric Total Variation Di- 
minishing (STVD) scheme of Yee [18]. In this scheme, gradients of characteristic 
variables are compared and selected such that no extraneous maxima or minima are 
introduced. This is accomplished through the use of a minmod function. This func- 
tion compares two values and returns the smallest in absolute magnitude if the signs 
are the same, or returns zero if the signs are different. The minmod portion of the 
second-order term in Eq. (3.5) is 

s mm = mmmod(si,S 2 ,s 3 ) . (3.13) 


The subscript 2 references the face at which the flux is being computed while 1 is the 
face behind and 3 is the face ahead. 

The implicit algorithm is therefore written as 
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.14) 


Here, e, f, and g are the normal fluxes in the £, ij, and ^-coordinate directions, re- 
spectively, while A, B, and C are their corresponding inviscid flux Jacobians. The 
absolute value Jacobians are computed using the Roe’s averaged variables and be- 
cause q is updated in planes parallel to the body, the scheme is point Jacobi within 
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the plane. This is a direct result of the fact that values used within the plane are 
all at the same time level. Values located above or below the plane being updated 
may be either the n or n+1 time level. These time levels are represented by the 
asterisk which remains in the formulation. For more details regarding the algorithm 
and relaxation strategies, one should consult Ref. 6. 


3.3 Boundary Conditions 

Figure 3.2 gives a schematic of a typical wall boundary. In this figure, (a) 
denotes the plane of cell-centered values a half cell above the body surface, (b) denotes 
the plane of cell-centered values a half cell below the body surface, and x denotes the 
location of points on the surface. The values at (b) are required for the computation 
of the first and second-order dissipation associated with the Roe and Yee methods, 
respectively. The values used at (b) are determined by simply equating them with 
the values used at (a). 

Values on the wall are determined such that surface tangency is observed. 
This method extrapolates values to the wall, and a wave correction is then performed 
on the values so that surface tangency is satisfied. In general, just extrapolating 
values to the surface will not meet this requirement. In order to determine the 
corrected values (c) at the wall, predicted values at the wall are found by first-order 
extrapolation using the values at (a). These values, along with the wave correction 
equations [19] 

Pa Pc 

pl pl 

U.- 2 -% = V,- 2-^-r , 

7-1 7-1 

u e = 0 , 

K = v c , 
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W a = W e , (3.15) 

are then used to determine the appropriate values at the wall. Here, U is the con- 
travariant velocity normal to the body surface, V and W are the contravariant veloc- 
ities tangent to the surface, and C is the speed of sound. 


3.4 Aerodynamic Quantities 


There are a number of different aerodynamic quantities which are used to 
describe the aerodynamic characteristics of a vehicle. Many of these quantities are 
based on the surface forces acting on the vehicle. The total surface force vector, F is 
given in integral form as 

F — J J P n da , (3.16) 

where P is the pressure on the surface, n is a unit vector normal to the surface, cr 
is an element of surface area, and S is the total surface area. The yaw, lift (L), and 
drag (D) are found by taking the x, y, and z components of F , respectively. The 
coefficient of lift, C/, and the coefficient of drag, C <j, are defined by 


Ci = 


= 


L/A 

\p« olKol 2 ’ 

D/A 

\poo\vL \ 2 ’ 


(3.17) 

(3.18) 


where v£> is the free stream velocity vector and A is the base area of the vehicle. The 
pitching-moment, M, for the vehicle is given by 


M = J J P (n x m) da , (3.19) 

where m is a vector extending from the vehicle’s center of gravity to the center of an 
element of area on the surface. The coefficient of pitching-moment is given by 

M/A 


C m — 


\Poo I H. 1 2 / 


(3.20) 
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where l is the length of the model. The coefficient of pressure, C p , is defined by 


Q, 


P-Pr, 


1 „ 

2 Poo 


V^ool 2 


( 3 . 21 ) 


where P ^ is the free-stream pressure. 



Chapter 4 


FLOWFIELD COMPUTATIONS AND 
DISCUSSION OF RESULTS 

4.1 Introduction 

Flowfield computations are performed on the (101x51x35) singleblock vol- 
ume grid using the LAURA code [6,8,1 1] described earlier. The Cray-2 computer 
at NASA Langley Research Center is used to make these computations. The code 
requires 13 Mwds of memory for a grid of 180285 points. Each iteration requires 
3.89 CPU seconds if the Jacobians are only updated every 20th iteration. Therefore, 
21.6 CPU microseconds are required per iteration per grid point. Computations are 
performed at Mach 6.0 and angles of attack of 0, 6, and 12 degrees. Free-stream 
conditions which existed in the 1968 study [1] are used in the LAURA code in order 
to make an accurate comparison between experimental and computational results. 

The base diameter (d in Figs. 2.10 and 2.12) of the 11.3° inner cone is 
assumed to be 4.0 inches. This creates a vehicle 8.0 inches in length with a base 
area of 14.85 square inches. These values are used as the characteristic length and 
characteristic area of the model, respectively (as was done in the 1968 study). The 
ratio of specific heats, 7, is taken to be 1.4 for these computations. 

Free-stream velocities of Mach 6.0 fall within the hypersonic range of fluid 
flow. A great deal of research has been done in the field of hypersonic flow over the 
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last few years. Blunt-nosed vehicles in particular have shown a number of common 
aerodynamic traits [12]. One of these is the formation of a strong normal shock in 
front of the vehicle. Immediately behind this shock is a region of subsonic flow with 
hypersonic and supersonic flow in the surrounding regions (Fig. 4.1). This type of flow 
is also characterized by a highly compressed region in front of the nose, overexpansion 
around the nose, and recompression downstream of the nose. Hypersonic flow about 
a cone is also well documented [12]. It is characterized by the formation of a strong 
shock wave coming off the cone with supersonic and hypersonic flow behind the shock. 

Because this study makes use of the Euler equations to determine the flow 
characteristics around the vehicle, diffusion and thermal conductivity effects are not 
accounted for. Values of drag and pitching moment should be below those deter- 
mined experimentally since only form drag, and not viscous drag, is computed by the 
Euler equations. Also, the change in enthalpy normal to the surface is assumed to be 
constant. Energy interaction by means of chemical reaction, radiation, molecular ro- 
tation, and molecular vibration are not accounted for within the flowfield. Therefore, 
temperatures on the surface found computationally are expected to be higher than 
those found experimentally. 


4.2 Cases Studied 

Euler computations are performed on the vehicle at three angles of attack 
(a = 0, 6, and 12 degrees) and assuming a free-stream Mach number of 6.0. 

Case 1: 0° angle of attack 

The contour plots of the coefficient of pressure (C p ) on the surface, symme- 
try planes, and exit plane are given in Fig. 4.2 for q = 0°. A strong normal shock wave 
forms just in front of the vehicle. The normal shock becomes an oblique shock and 
continues to the exit plane where the conical shape of the shock becomes apparent. 
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Fig. 4.1 Hypersonic How around a blunt body. 
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The pressure contours throughout the flowfield are symmetrical, which is expected for 
an angle of attack of 0°. It is interesting to note a region of high pressure surrounded 
by regions of lower pressure on the outer cones of the forebody. An enlarged view 
of the C p contours on the first thirty tj lines of the surface is also given in Fig. 4.2. 
The black region of the nose indicates the very high pressure gradient which is ex- 
pected in this region. Also demonstrated in this figure are the low pressure regions 
slightly downstream of the nose. They result from the overexpansion out of the high 
pressure region in front of the nose. The coefficients of lift, drag, and pitching mo- 
ment found computationally are -2.776x10-*, 0.1373, and 4.128xl0- 5 , respectively, 
while these values found experimentally are 0.0, 0.145, and -5.0x10 4 , respectively. 
This computed lift is near 0.0 computationally which is expected for 0° angle of at- 
tack. The magnitude of the coefficients of drag and pitching moment, however, are 
underestimated because viscous effects are not accounted for. 

Case 2: 6° angle of attack 

A plot of C p contours is given in Fig. 4.3 for a = 6°. This figure shows a 
strong shock on the windward side of the vehicle. A close examination of Fig. 4.3 
compared to Fig. 4.2 shows that the windward shock has moved slightly closer to the 
body while the leeward shock has moved slightly away. An enlarged view of the nose 
of the vehicle shows closed pressure contours on the windward side. These contours 
are again the result of overexpansion near the nose. The coefficients of lift, drag, 
and pitching moment found computationally are 0.2455, 0.1659, and -1.001x10 , 

respectively, while these values found experimentally are 0.245, 0.18, and -2.5x10 , 

respectively. The coefficient of lift compares very well with the experimental values 
but the coefficient of drag and pitching moment are below those found experimentally, 


as expected. 
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Case 3: 12° angle of attack 

A plot of C p contours is given in Fig. 4.4 for a = 12°. A strong shock 
has formed on the windward side of the vehicle. When compared to the a = 6° 
case, the windward shock has moved closer to the vehicle while the leeward shock 
has moved farther away from the vehicle. There are also large pressure gradients 
located on the outer cones of the forebody. These result from the fact that the outer 
cones are becoming leading edges at higher angles of attack. The forebody /afterbody 
junction seems to cause large pressure gradients at higher angles of attack as well. The 
enlarged view of the model shows the closed contours of low pressure near the nose 
which are caused by overexpansion in this region. The low pressure has continued 
to move towards the nose on the windward side as the angle of attack has increased. 
The coefficients of lift, drag, and pitching moment found computationally are 0.4755, 
0.2727, and -4.739xl0 -3 , respectively, while these values found experimentally are 
0.482, 0.30, and -5.0xl0 -3 , respectively. The coefficient of lift for this case and the 
two other cases, compares very well with the experimental results. The coefficient of 
drag and pitching moment, however, are again below those found experimentally. 

4.3 Comparisons 

Figures 4.5 through 4.8 illustrate the computational and experimental values 
of the coefficients of lift, drag, and pitching moment, as well as lift -over-drag for the 
three cases studied. Figure 4.5 in particular shows that values of the coefficient of lift 
found computationally compare very well with those found experimentally. Figure 4.6, 
however, shows that the magnitudes of the coefficient of drag found computationally 
are less than those found experimentally. This is due to the fact that viscous effects are 
not taken into account. The magnitudes of the coefficient of pitching-moment about 
the x-axis found computationally are also less than the values found expermentally. 
This is again the result of viscous effects not being taken into account. Corrected 
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Fig. 4.4 C p contours for Mach 6.0, a = 12°. 
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values for the coefficients of drag and pitching-moment are plotted in Figs. 4.6 and 
4.7, respectively. They are found by adding the difference between the computational 
and experimental values at 0° angle of attack to the remaining computational values. 
This gives a slight correction to allow for viscous effects. Figure 4.8 compares the 
computational lift-over-drag ratios to those found experimentally. Since drag is un- 
derestimated computationally, lift-over-drag is overestimated. The lift-over-drag at 
0° angle of attack, however, is accurate since the value of lift is virtually zero. 

Figure 4.9 gives a plot of the density along the stagnation line of the vehicle 
at 0° angle of attack. This figure shows a strong normal shock coming off the vehicle. 
The standoff distance of this shock is approximately 0.03 inches. Figures 4.10, 4.11, 
and 4.12 give the distribution of C p along the top, bottom, and side of the model, 
respectively. Figure 4.10 shows that as the angle of attack increases, the pressure 
on the top of the vehicle (which is on the leeward side at positive angles of attack) 
decreases. At 0° angle of attack, a slight dip in the pressure distribution is seen near 
the nose. This is a result of the overexpansion which occurs in this region. A reduction 
in pressure is also seen at 4.0 inches which is the location of the forebody/afterbody 
junction. This is the result of an expansion which occurs in this region. Figure 
4.11 gives a plot of C p along the bottom of the vehicle. This figure illustrates that 
as the angle of attack increases, the pressure increases. The overexpansion near 
the nose becomes greater on the bottom of the vehicle at higher angles of attack. 
The expansion at the forebody /afterbody junction is also greater at higher angles of 
attack. Figure 4.12 shows the variation of C p along the side of the vehicle for 0° angle 
of attack. The overexpansion and recompression near the nose is clearly seen in this 
figure. Also, since there is a greater discontinuity in slope on the side of the model, 
a greater expansion occurs there. This expansion is shown by the large reduction in 
pressure at this point. 











Chapter 5 

Conclusion 


The purpose of this study has been to computationally verify the results 
obtained in 1968 and to continue the study of this vehicle. This study has also, 
however, helped to further verify the LAURA code. The results obtained agree very 
well with the experimental results. The lift calculated by the LAURA code agrees 
almost entirely with the lift computed experimentally. The drag and pitching moment 
are slightly below the experimental results but they follow the expected trends. The 
vehicle has a maximum lift-over-drag ratio of 1.75 computationally as compared to a 
value of 1.57 found experimentally. The pressure distribution on the surface shows 
that there is a highly compressed region in front of the nose, overexpansion around 
the nose, and recompression downstream of the nose. An expansion also occurs at 
the forebody/afterbody junction where there is a change in slope. This expansion is 
pronounced on the side of the vehicle because of the greater change in slope. 

There are currently plans for the continuation of the research on this vehicle. 
The effects of different angles of attack and Mach numbers will be examined. Also, 
new vehicles will be constructed with higher degrees of nose bluntness. These vehicles 
will be tested using different methods to determine the effects that nose bluntness 
has on the flow characteristics of the vehicle. 
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PROGRAM ZFIND 

************************************************************** 

* JOHN STEWART FEBRUARY 10, 1990 * 

* PROGRAM WHICH FINDS THE Z-VALUES FOR EACH CROSS-SECTION * 

* TO BE USED IN PROGRAM SHIP. CALCULATIONS ARE BASED * 

* ON ARCLENGTH OF INNER CONES AND SPHERICAL NOSE. * 

* * 

* NI = NUMBER OF CROSS-SECTIONS TO BE FOUND (MUST * 

* MATCH NI IN PROGRAM SHIP) * 

* NS = NUMBER OF ARC LENGTH DIVISIONS USED FOR CALC. 'S* 

+ ZI = LENGTH OF AFTERBODY * 

* ZII = LENGTH FROM INNER CONE TANGENTS TO AFTERBODY * 

* 21 1 1 = LENTH OF SPHERICAL NOSE * 

* D = MAXIMUM AFTERBODY INNER CONE DIAMETER * 

* STMAX = MAXIMUM ARC LENGTH * 

************************************************************** 
PARAMETER (NI=200,NS=5000) 

DIMENSION Zl(NS) ,AL(NS) , Z(NI+1), Y(NI+1) 

REAL KSP 

OPEN ( 6 , F I LE= ' ZDATA ' , FORM= ' FORMATTED ' ) 

D=4 . 0 
ZI=D 

ZI I=. 97039*D 
ZII I=D-ZI I 

DEFINE THREE POINTS FOR THREE EXPONENTIAL CURVES TO FIT 
SMOOTHLY THROUGH. POINTS ARE (TAU*,AA*) . DISTRIBUTIONS 
OF ARC LENGTH WILL BE CHOSEN OFF AA AXIS FOR EQUAL DIVISIONS 
OF ARC LENGTH ON TAU AXIS. KSP CONTROLS THE AMOUNT OF 
CURVITURE OF THE CURVES. 

AA1=. 02519 
AA2= . 25 
AA3=. 509 
TAU1=. 13 
TAU2=. 35 
TAU3=. 7 
KSP= . 01 

STMAX = 2 . 07703*D 

CALL SPACE ( AA1 , AA2 , AA3 , TAU1 , TAU2 , TAU 3 , KSP . NI ♦ 1 . STMAX , Y ) 

LOOP FINDS THE VALUES OF ARC LENGTH FOR EACH OF NS 
CROSS-SECTIONS BASED ON DISTRIBUTION FOUND ABOVE 

DO 31 1=1, NS 

21 ( I > = 2 . 0*D* ( FLOAT ( 1-1 )/FLOAT(NS-5 ) ) 

IF (Zl(I).LT. .02961*D) THEN 
AL( I ) = . 04*D*AC0S( ( . 04*D-Z1 ( I ) )/( . 04*D) ) 

GO TO 31 

ELSE IF (Zl(I).LT. (1.0*D)) THEN 
AL( I ) = . 05231*D+1 . 035616*(Z1(I)-ZIII) 

GO TO 31 
ELSE 

AL( I )=1.05726*D+(Z1( I ) -ZI I-ZI 1 1 ) *1 . 019769 
GO TO 31 
END IF 
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31 CONTINUE 


FINDS CROSS-SECTIONAL Z LOCATION BASED ON ARC LENGTH 
DISTRIBUTION 


K=1 

DO 26 I = 1, NI+1 
28 k = K ♦ 1 

IF(Y( I ) .GT. AL(K) ) THEN 
GO TO 28 
ELSE 

Z(I) =( ( Y( I ) -AL( K- 1) )*(Z1(K)-Z1(K-1) ) )/( AL(K) -AL(K- 1 ) ) + Z1(K 
K = K -1 
END IF 

IF(I.GT.l) WRITE ( * , * ) I, Z(I)-Z(I-1).Z(I) 

26 CONTINUE 

WRITES OUT RESULTS TO A FILE ZDATA WHIRE I=NI+1 IS THE 
NOSE OF THE SPACECRAFT 

DO 1 I = NI+1 , 1. -1 
WRITE ( 6 , * ) Z(I) 

1 CONTINUE 
STOP 
END 

SUBROUT INE SPACE(A1,A2,A3,TAU1, TAU2 , TAU3 . K2 , N , ST . S 1 ) 
************************************************************** 

* SUBROUTINE WHICH FITS A SMOOTH CURVE THROUGH THREE POINTS 

* USING THREE EXPONENTIAL EQUATIONS OF THE FORM: 

* X = ( EXP ( KSP*X ) - 1 ) / (EXP(KSP)-l) 
************************************************************** 

DIMENSION SI (N) 

REAL Kl.K2,K3,K4,K5,NOM, IX 
I COUNT=0 
K3=-K2 

300 CONTINUE 

TEMP=A1*K2* (TAU2-TAU1 )/( ( A2-A1 ) *TAU1 ) 

DETAL=1. /( 1-EXP ( -K2 ) ) 

DETAR=K3/ ( EXP ( K3 ) - 1 . ) 

C=TEMP * DETAL 
DF=DETAR-C 

DFDK=(EXP(K3)-1.-K3*EXP(K3) )/( EXP ( K3 ) -1 ) **2 
DK3=-DF/DFDK 

IF(ABS(DK3) .LT. .00001) GO TO 200 
K3=K3*DK3 
I COUNT= I COUNT* 1 
IF( ICOUNT.GT. 20) GO TO 200 
GO TO 300 
200 CONTINUE 
K4=-K3 

301 CONTINUE 

TEMP=A2 * K3 * ( TAU3 - TAU2 ) / ( ( A3 - A2 ) * TAU2 ) 

DETAL=1. /( 1-EXP ( -K3 ) ) 

DETAR=K4/ ( EXP ( K4 ) - 1 . ) 

C=TEMP*DETAL 



DF=DETAR-C 

DFDK= ( EXP ( K4 ) - 1 . - K4* EXP ( K4 ) )/( EXP( K4 ) - 1 ) **2 
DK4=-DF/DFDK 

IF(ABS(DK4) . LT. . 00001 ) GO TO 201 

K4=K4+DK4 

GO TO 301 

201 CONTINUE 
K5=-K4 

302 CONTINUE 

TEMP=A3*K4* ( 1 . -TAU3 )/( ( 1 . -A3 ) *TAU3 ) 

DETAL=1 . /( 1-EXP ( -K4 ) ) 

DETAR=K5/ ( EXP ( K5 ) - 1 . ) 

C=TEMP *DETAL 
DF=DETAR-C 

DFDK= ( EXP ( K5 ) - 1 . - K5 * EXP (K5) )/(EXP(K5)-l. )**2 
DK5=-DF/DFDK 

IF(ABS(DK5) .LT. .00001) GO TO 202 
K5=K5 +DK5 
GO TO 302 

202 CONTINUE 

1 CONTINUE 
SCALEX=ST 
SCALEF=SCALEX 
DO 10 1=1, N 

ETA= FLOAT ( I - 1 ) /FLOAT ( N- 1 ) 

IF(ETA.GE.TAUl) GO TO 2 

TERM=K2/TAU1 

NOM=EXP ( TERM* ETA ) - 1 . 

DNOM=EXP ( K2 ) - 1 . 

F=A 1 * ( NOM/DNOM ) 

GO TO 3 

2 CONTINUE 

I F ( ETA . GT . TAU2 ) GO TO 4 
TERM=K3/ (TAU2-TAU1 ) 

N0M=EXP(TERM*ETA-TERM*TAU1) -1. 

DNOM=EXP ( K3 ) - 1 . 

F=A1 + ( A2 - A1 ) * NOM/DNOM 
GO TO 3 

4 CONTINUE 

I F ( ETA . GT . TAU3 ) CO TO 5 
TERM=K4/ ( TAU3 -TAU2 ) 

NOM=EXP ( TERM*ETA-TERM*TAU2 ) - 1 . 

DNOM=EXP ( K4 ) - 1 . 

F=A2 ♦ ( A3 - A2 ) * NOM/DNOM 
GO TO 3 

5 CONTINUE 
TERM=K5/( 1 . -TAU3 ) 

NOM=EXP ( TERM* ETA-TERM*TAU3 ) -1 . 

DNOM=EXP ( K5 ) - 1 . 

F=A3* ( 1-A3 ) * NOM/DNOM 

3 CONTINUE 
IX=SCALEX*ETA 
Sl( I)=SCALEF*F 

IF ( I . EQ. ( 35+1 )/2 ) WRITE( * , * ) ' ***** ' , IX. Sl( I ) 
IF (I.EQ. (37+l)/2) WRITE( * , * ) '******'. IX. Sl( I ) 
10 CONTINUE 
RETURN 
END 
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PROGRAM SHIP 


★ 

★ 

★ 

★ 

★ 

* NOTE: 

* 


JOHN STEWART FEBRUARY 10, 1990 

PROGRAM TO ANALYTICALLY DEFINE THE SURFACE GRID 
OF THE SPACESHIP FOUND IN NASA TECHNICAL NOTE: 

NASA TN D- 4598 , JUNE, 1968. 

IN PROGRAM, SHIP IS ORIENTED SUCH THAT THE NOSE IS 
AT (0,0,0) AND THE 2- AXES RUNS DOWN THE LENGTH OF 
THE SHIP. THE X-AXIS INTERSECTS THE CYLINDER AND 
OUTER CONE CENTERS. 1=1 IS THE BACK OF THE SHIP 
AND I=NI ♦ 1 IS THE NOSE OF THE SHIP. 

OUTPUT DATA IS ORIENTED SUCH THAT X=Z, Y=X, Z=Y AND 
1=1 IS AT THE NOSE. 


* 

* 

* GIVEN DATA: 

* 

★ 

* 

★ 

★ 

* 

★ 

* 

* 

★ 

★ 

* 

* 

* 

* 

★ 

★ 

* 

* 

* 

* 

* 


CONANGl = AFTERBODY CONE ANGLE (RADIANS) 

CONANG2 = FOREBODY CONE ANGLE (RADIANS) 

D = MAXIMUM DIAMETER OF AFTERBODY CONE 

ELLD = DISPLACEMENT BETWEEN CONE EDGE AND CYLINDER CENTER 
IN AFTERBODY MEASURED PERPENDICULAR TO CONE WALLS 
ZI = LENGTH OF AFTERBODY 

ZII= LENGTH FROM INNER CONE TANGENTS TO AFTERBODY 
ZIII= LENCTH FROM SPHERE TO INNER CONE TANGENTS 
ZIV= LENGTH OF SPHERICAL NOSE 
RPIPE = RADIUS OF PIPES IN AFTERBODY MEASURED 
PERPENDICULAR TO CONE WALLS 

RDIV = RADIUS OF CONE AT FOREBODY-AFTERBODY CONNECTION 
ELLD1 = ELLIPSE DISPLACEMENT ON FOREBODY AT CONNECTION 
B1 = MAJOR AXIS OF ELLIPSE ON FOREBODY AT CONNECTION 
A1 = MINOR AXIS OF ELLIPSE ON FOREBODY AT CONNECTION 
XO = RADIUS OF NOSE IN XY PLANE AT INNER NOSE TANGENTS 
XOl = RADIUS OF NOSE IN XY PLANE AT OUTER CONE TANGENTS 
NI = NUMBER OF DIVISIONS DOWN LENGTH OF SHIP 
N = NUMBER OF POINTS OUTLINING SHIP IN XY-PLANE PER QUAD. 
L - NUMBER OF GRID POINTS TO LEAVE OFF THE SURFACE GRID 
AT THE NOSE. 


************** **************************************************** 
PARAMETER (NI =200. N=50) 

PARAMETER ( NN=N+ 1 , L=0 ) 

DIMENSION Z( 1+NI-L) , A(NI+1-L) 

♦,B(NI+1-L) , R(NI+1-L) 

♦,XN(NI+1-L), CN(NI+1-L) , ED(NI+1-L), CC(NI+1-L) 

♦ , XC (NI+l-L), Y(NI+1-L,NN*2-1) 

+.X(NI-L+1,NN*2-1) , XG (NI+l-L) 

♦, YN(NI+1-L) , DN(NI+1-L), R2(NN) 

+,B5(NI+1-L) , ANG(NN) 

REAL M(NI+1-L) , KSP. Kl, K2 

OUTPUT FILE FORMATTED FOR USE WITH PLOT3D ON THE IRIS 
WORKSTATION 


OPEN( 6, FILE= ' SHIPO' , FORM= ' FORMATTED ' ) 

INPUT FILE GIVING THE CROSS-SECTIONAL Z LOCATIONS 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


STARTING FROM THE BACK OF THE AFTERBODY AND MOVING FORWARD 
OPEN ( 7 , F I LE= ' 2DATA ' , FORM= ' FORMATTED ' ) 

OUTPUT FILE FORMATTED AS X, Y, 2 WHRE X IS NOW THE AXIS 
OF THE SPACESHIP AND N VALUES OF Y AND Z ARE GIVEN FOR 
EACH VALUE OF X 

OPEN ( 9 , FILE= ' SHI POD ' , FORM= ' FORMATTED ' ) 

OUTPUT FILE GIVING THE Y-LOCATION OF THE TANCENCY POINTS 
ON THE INNER CONE AND OUTER CONE ( FOREBODY ) OR 
CYLINDER (AFTERBODY) . 

OPEN ( 10 , FI LE= ' XNCN ' , FORM= ' FORMATTED ' ) 

C0NANG1 = .19722 
C0NANG2 = .26302 
D = 4.0 

XO = . O30624*D 
X01 = . 036072*D 
XODI FF=X0-X01 
ZI = D 

ZII = . 9704*D 

ZIII = ( .0296-.022723)*D 

ZIV = . 022723 *D 

RPIPE = . 12 *D 

RDIV = . 30018*D 

ELLD = . 08*D 

ELLDISP = ELLD/COS ( CONANG 1 ) 

ELLD1 = . 07146*D 
B1 = . 13249*D 
A1 = . 1203*D 

READ IN THE CROSS-SECTIONAL AXIAL LOCATIONS BEGINNING 
FROM THE BACK OF THE SPACESHIP 

DO 75 I = 1* NI+l-L 
READ( 7 , * ) Z ( I ) 

75 CONTINUE 
11=0 

CALCULATIONS OF CONE RADIUS (R) # MAJOR (B) AND 
MINOR (A) AXIS OF ELLIPSES CREATED BY PIPES, AND DISTANCE 
TO CENTER OF ELLIPSE (XC) FOR AFTERBODY. XG AND CG ARE 
THE POINTS OF TANGENCY GUESSES FOR SUBROUTINE POINTS. 

76 11=11+1 

IF (Z( II) .LT.ZII+ZIII+ZIV) GO TO 79 
ED (II) = ELLDISP 
A( I I ) = RPIPE 
B ( I I ) = A( II ) /COS ( CONANG 1) 

R( 1 1 ) = D/2 . + ( Z( 1 1 ) - ZI-ZII-ZIII-ZIV) * TAN ( CONANG 1 ) 

XC (II) = R( I I ) + ED (II) 

CG(II) = XC (II) + . 5*B( I I ) 

Z1=(CG(II)-XC(II))**2 
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Z2=Zl/( B( 1 1 ) *B( 1 1 ) ) 

Z3=(A( II )*A( II ) )/(B( II )*B( II )*B( II )*B( II) ) 

Xl=( Z1*R( II )*R( II )*Z3 )/( 1.-Z2) 

X2=( 1 . ♦ (Xl/( R( I I ) *R( I I ) ) ) ) 

XG( II ) = SQRT ( X1/X2 ) 

GO TO 76 

SAME CALCULATIONS FOR FOREBODY 
79 1111=11-1 

77 B ( 1 1 ) = ( Bl* ( Z ( II )-ZIV) )/(ZII+ZIII ) 

ED(II) = ELLD1* (Z( II )-ZIV)/(ZII+ZIII ) 

R( I I ) = ( ( (RDIV)-XO)*(Z( II )-ZIII-ZIV) )/ZII + XO 
XC(II) = ( ( ( RDIV+ELLD1 ) -XOl )*(Z( II )-ZIV) )/(ZII+ZIII ) +X01 
A( 1 1 ) = (A1*(Z( II )-ZIV) )/(ZII+ZIII ) 

RADIUS FOUND IN REGION ZIII WHERE Z IS NOW ON THE 
SPHERICAL NOSE 

IF(Z( II ) . LT. 2III+ZIV) THEN 
S5=. 04*D - Z ( 1 1 ) 

R( II )=SQRT( . 0016*D*D-S5*S5 ) 

ENDIF 

CG(II) = XC(II) + . 5*B( II) 

Z1=(CG( 1 1 ) -XC( II) ) **2 
Z2=Zl/( B( 1 1 ) *B( 1 1 ) ) 

Z3=(A(II)*A(II))/(B(II)+B(II)*B(II)*B(II)) 

X1=(Z1*R( I I )*R( I I ) *Z3 )/( 1 . -Z2 ) 

X2=( l.+(Xl/(R(II)*R(II)))) 

XG (II) = SQRT(X1/X2 ) 

11=11+1 

IF (II.GT. (NI+l-L)) GO TO 2 
IF (Z(II) .LT.ZIV) GO TO 2 
IF ( XC (II-1)+B(II-1). LT .R(II-l)) GO TO 2 
GO TO 77 
2 CONTINUE 
III = 1 1 -1 
XX1=1 .0 

WRITE( 10, * ) I XN J AFTER XN CN J AFTER CN 

&YMAX ' 

P111=0 THE MAXIMUM OF 15 ITERATIONS IN SUBROUTINE POINTS 
HAS NOT BEEN ACHIEVED ( NEWTON -RAPHSON) 

Plll=l THE MAXIMUM HAS BEEN ACHIEVED SO THE PROGRAM USES 
SUBROUTINE POINTSB FROM NOW ON BECAUSE IT IS MORE 
EFFICIENT (BISECTION). 

P111=0.0 

DO 3 I = 1, III 

CALCULATING XN AND CN -- TANGENCY POINTS 
IF (P111.EQ.0.) THEN 

CALL POINTS(XC( I),A(I),B(I),R(I), XG( I ) , CG( I ) , Pill , 

+XN(I).CN(I)) 

ENDIF 
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IF (Plll.EQ. 1. ) THEN 

CALL POINTSB ( XC ( I ) , A ( I ) f B ( I ) , R { I ) , XN ( I ) , CN ( I ) ) 

WRITE ( * » * ) '****' 

ENDIF 

WRITE ( * , * ) I, XN ( I ) , CN ( I ) 

C 

C EQUATIONS TO CALCULATE THE ACTUAL POINTS WHICH WILL 

C OUTLINE THE SHIP IN THE XY PLANE IN THE FIRST QUADRANT. 

C YN IS THE Y- LOCATION OF TANGENCY POINT XN 

C DN IS THE Y- LOCATION OF TANGENCY POINT CN 

C M IS THE SLOPE OF THE LINE JOINING THE POINTS 

C B5 IS THE Y INTERCEPT OF THE LINE 

C 

YN{ I ) = SQRT(R(I)*R(I)-XN(I)*XN(I)) 

DN ( I ) = A ( I ) * SQRT ( 1 . 0 - ( ( ( CN ( I ) - XC ( I ) ) /B (I))**2)) 

M( I ) = ( YN ( I ) - DN (!))/( XN ( I ) - CN ( I ) ) 

B5 ( I ) = YN(I) - M ( I ) * XN ( I ) 

3 CONTINUE 

X AND Y POINTS ARE CHOSEN AT CONSTANT ANGLES FOR EVERY J. 
SUBROUTINE SPACE2 USES TWO EXPONENTIAL CURVES WHICH 
CONNECT AND ARE Cl CONTINUOUS AT TAU. EQUATION OF CURVES IS: 
X= ( EXP ( X*K1 ) - 1 ) / (EXP(Kl)-l) 

Kl=. 5 
K2=2 . 0 
TAU- . 5 

ANGMAX=1. 5707963 

CALL SPACE2 ( K1 , K2 , TAU, NN, ANGMAX , ANG) 

DO 37 J = 1, NN 
ANG(J) = 1 . 5707963 -ANG ( J ) 

37 CONTINUE 

SETTING THE X AND Y VALUES TO THEIR EXACT VALUES 

DO 29 I » 1,111 
X(I,1)=0 
Y( I , 1 ) =R{ I ) 

X(I,NN)=XC(I)+B(I) 

Y( I , NN)=0 . 0 

BISECTION METHOD TO FIND THE VALUES OF X AND Y WHICH 
CORRESPOND TO THE SPECIFIED ANGLES 

DO 40 J - 2 , N 
XS1=0.0 
XS2=XC ( I ) +B( I ) 

N5“0 

30 XG1= ( XS2 +XS 1 ) /2 . 

N5=N5+1 

CALL YFIND(XGl t XN(I),CN(I) B R( I ) f A( I ) ,B( I ) ,XC( I ) ,M( I ) , B5( I ) .YG1) 
TEST=ATAN ( YG1/XG1 ) 

REP=TEST- ANG ( J ) 

IF(REP.LT.O) THEN 
XS2=XG1 
ELSE 
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XS1=XC1 

ENDIF 

I F ( N5 . GT . 20 ) GO TO 41 
GO TO 30 
41 X( I , J) =XG1 
Y( I , J ) =YG1 
40 CONTINUE 
29 CONTINUE 

DEFINING THE NOSE OF THE SPACESHIP TO BE POINT (0,0,0) 
AND DEFING THE SPHERICAL SHAPE OF THE NOSE (REGION 2IV) 

78 DO 81 J = 1, NN 

IF (Z(II).EQ.O.O) THEN 
Z( II )=0. 

X(II.J)= o. 

Y(II,J)= 0. 

GO TO 81 
ENDIF 

S5 = . 04*D - Z ( 1 1 ) 

R5 = SQRT( . 0016*D*D - S5*S5) 

X(II,J) = R5*COS( ANG( J) ) 

Y ( 1 1 , J ) = R5* SIN( ANG( J ) ) 

81 CONTINUE 
11=11+1 

IF ( II .EQ.NI+2-L) GO TO 80 
GO TO 78 
80 CONTINUE 

ONLY ONE QUARTER OF THE MODEL WAS CONSTRUCTED SO THE 
POINTS ARE EXTENDED TO CREATE HALF A MODEL. 

13 DO 9 I = 1, NI+l-L 
DO 9 J = 1, N 
X ( I , NN* 2 - J ) =X ( 1 1 J ) 

Y( 1 1 NN*2- J )=-Y ( I * J ) 

9 CONTINUE 

PLOT3D FILE IS CREATED 

WRITE (6,*) NI +1-L, NN*2-1, 1 
DO 55 J * 1, NN*2-1 
DO 55 I = 1, NI+l-L 
WRITE ( 6 , * ) Z(I) 

55 CONTINUE 

WRITE( 6, * ) X 
WRITE ( 6 , * ) Y 

WRITE( 9 , * ) NN*2-1, NI-L+1 
DO 57 I = NI-L+1 , 1, -1 
DO 57 J = 1, NN*2- 1 
WRITE( 9 , * ) Z ( I ) f X( I , J ) , Y{ I , J) 

57 CONTINUE 
PPP=0 . 

WRITING OUT THE XNCN FILE TELLING WHERE THE LOCATIONS 
OF THE TANGENCY POINTS ARE FOR EACH CROSS-SECTION 
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DO 204 I = III, 1, -1 
DO 205 J = 1, 51 
I F( PPP . EQ . 2 . ) GO TO 205 
IF(PPP. EQ. 1 . ) GO TO 206 
I F ( X ( I , J ) . GT . XN ( I ) ) THEN 
PPP=1 . 

JJ1=J 
END IF 

206 I F ( X ( I , J ) . GT . CN ( I ) ) THEN 
PPP=2 . 

JJ2-J 
END IF 

205 CONTINUE 

WRITE ( 10 , 203 ) 202-1, XN( I ) , JJ1, CN(I), JJ2 , XC(I)+B(I) 

PPP=0. 

204 CONTINUE 
STOP 

200 FORMAT ( 5X, 12 , F10 . 3 , F10 . 3 , F10 . 3 ) 

203 FORMAT ( 5X , 13 , 5X, F9 . 6 , 5X. 12 , 5X, F9 . 6 , 5X, 12 . 5X t F9 . 6 ) 

END 

SUBROUTINE POINTS (XC , A , B , R ( XG, CG , PI 1 1 , XN, CN) 

ft************************************************************** 

* SUBROUTINE TO DETERMINE THE POINTS WHERE A LINE IS TANGENT* 

* TO BOTH A CIRCLE OF RADIUS R, AND AN ELLIPSE WITH MAJOR * 


* AXIS B AND MINOR AXIS A USING NEWTON- RAPHSON METHOD. * 

* XG = INITIAL GUESS FOR X ON CIRCLE * 

* CG = INITIAL GUESS FOR X ON ELLIPSE * 

* XN = X ON CIRCLE * 

* CN = X ON ELLIPSE * 


P1=0. 

C = CG 

X = XG 
P=l. 

J=0 

10 J=J+1 

IF (X.GT.R) X=R* . 999 
IF (C.LT.XC) C=XC*1.001 
IF(J.EQ.15) THEN 
Plll=l. 

RETURN 
END IF 

Z1 = (C - XC) 

Z2 = Zl/B 

Z3 = ( (R**2)-(X**2) ) 

F = (Z1/SQRT(1-(Z2**2) ) )*(A/(B**2) )-X/SQRT(Z3) 

FX = ((A/(B**2))*Z1)/((1-(Z2**2))**.5)-(X**2)/(Z3**1.5) 
+-1.0/(Z3**.5) 

FC = (A/(B**3))*(Z1**2)/((1.-(Z2**2))**1. 5 )-X/(Z3** . 5 ) + 
+(A/(B**2) )/((l-(Z2**2))**.5) 

C = 2.*(R**2)-2.*C*X-2*A*(((1.-(Z2**2))*Z3)**.5) 

GX = -2.*C+(2*X*A*(1.-(Z2**2)))/(((1.-(Z2**2))*Z3)**.5) 
GC a -2 . *X+ ( ( ( 2 . *A)/B) *Z3*Zl)/( ( (l.-(Z2**2) )*Z3)**.5) 

XI = X - ( F*GC-FC*G)/( FX*GC-FC*GX) 
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C5 - C - ( F*GX-FX*C)/( (FC*GX-FX*GC)*P) 

IF ( (C5-XC)/B.GE. 1. ) THEN 
P=100. 

GO TO 10 
ENDIF 

XG1 = ABS(Xl-X) 

CGI = ABS(C5-C) 

IF (XG1 . GT . ABS ( . 005*X1 ) ) THEN 
X=X1 
C=C5 

GO TO 10 

ELSE IF (CGI . GT. ABS ( . 005*C5 ) ) THEN 
X=X 1 
C=C5 

GO TO 10 
ENDIF 
XN = XI 
CN = C5 
RETURN 
END 

SUBROUTINE YFINDfX. XN. CN. R. A, B, XC . M. B5 . Y) 

* SUBROUTINE TO FIND Y CIVEN X * 

real*m*********************************************************** 

IF (X.LE.XN) THEN 
Y=SQRT(R*R-X*X) 

GO TO 25 
ENDIF 

IF(X.GE.CN) THEN 
Y = A*SQRT(1.0-( ((X-XC)/B)**2) ) 

GO TO 25 
ENDIF 

Y a M*X + B5 
25 RETURN 
END 


C 

C 

C 


SUBROUTINE XFIND(Y, YN,DN,R,A,B,XC,M,B1,X) 

: i»™«******** ;*;**;*;*i;?2********* ******** 

SUBROUTINE FINDS X CIVEN Y * 

- — — — *************** ******** 

IF (Y.LE.DN) THEN 
X = B*SQRT( 1 . 0- ( ( Y/A) **2 ) ) + XC 
GO TO 30 

ELSE IF (Y.GE.YN) THEN 
X = SQRT( R*R - Y*Y) 

GO TO 30 
ELSE 

X = ( Y-Bl )/M 
ENDIF 
30 RETURN 
END 


SUBROUTINE POINTSB(XC, A, B. R, XN. CN) 
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c *************************************************************** 

C * SUBROUTINE TO DETERMINE THE POINTS WHERE A LINE IS TANGENT* 

C * TO BOTH A CIRCLE OF RADIUS R, AND AN ELLIPSE WITH MAJOR * 

C * AXIS B AND MINOR AXIS A BY USING THE BISECTION METHOD * 

C * Cl = INITIAL GUESS FOR X * 

C * C2 = INITIAL GUESS FOR X * 

C * XN = X ON CIRCLE + 

C * CN = X ON ELLIPSE * 

c *************************************************************** 

C1=XC 
C2=XC+B 
J= 0 

Z3=(A*A)/(B*B*B*B) 

10 J=J+1 

C = ( C1+C2 )/2 . 

Z1=(C-XC)**2 
Z2=Zl/( B*B ) 

X1=(Z1*R*R*Z3)/(1. -Z2) 

X2=( 1. +(X1/(R*R) ) ) 

X=SQRT { X1/X2 ) 

G = R*R - C*X - A*SQRT( ( 1 . -Z2 ) * (R*R-X*X) ) 

G = ( ( C-XC ) *A )/( ( SQRT ( 1-Z3 ) ) *B*B ) - ( X/( SQRT( Z1 ) ) ) 

IF (G.GT.0.0) THEN 
C1=C 
ELSE 
C2=C 
END IF 

IF ( J. LT.20) GO TO 10 
XN = X 
CN = C 
RETURN 
END 

SUBROUT I NE SPACE2 ( K1 , K2 , TAU , N , ANGMAX , ANG ) 

***************** ********************************************* 

* SUBROUTINE USES TWO EXPONENTIAL CURVES TO GIVE A * 

* DISTRIBUTION FOR THE ANGLES USED IN PROGRAM * 

DIMENSION ANG(N) 

REAL K1 , K2 , K3 , NOM 
I COUNTED 
K3=-K2 
300 CONTINUE 

TEMP=K1*K2* ( l-TAU)/( ( 1-K1 ) *TAU) 

DETAL=1 . /( 1-EXP ( -K2 ) ) 

DETAR=K3/ ( EXP ( K3 ) - 1 . ) 

C=TEMP * DETAL 
DF=DETAR-C 

DFDK=( EXP(K3 ) -1 . -K3*EXP(K3 ) )/(EXP(K3 ) -1 ) **2 
DK3=-DF/DFDK 

IF(ABS(DK3) .LT. .00001) GO TO 200 
K3=K3+DK3 
I COUNT= I COUNT* 1 
I F ( ICOUNT.GT.20) GO TO 200 
GO TO 300 
200 CONTINUE 


SCALEX=ANCMAX 

SCALEF=SCALEX 

1 CONTINUE 

DO 10 1=1, N 

ETA=FLOAT ( I - 1 ) /FLOAT ( N- 1 ) 
IF(ETA.CE.TAU) GO TO 2 
TERM=K2/TAU 
NOM=EXP ( TERM*ETA ) - 1 . 

DNOM=EXP ( K2 ) - 1 . 

F=K1 * ( NOM/DNOM ) 

GO TO 3 

2 CONTINUE 
TERM=K3/( 1-TAU) 

NOM=EXP ( TERM* ETA-TERM*TAU ) -1. 

DNOM=EXP ( K3 ) - 1 

F=K1 + ( 1 , -K1 ) * NOM/DNOM 

3 CONTINUE 
IX=SCALEX*ETA 
ANG ( I ) =SCALEF*F 

10 CONTINUE 

11 CONTINUE 
RETURN 
END 




